Generalized traveling wave in heterogeneous 

media: 

Allee effect can inhibit invasion 



benoit . sar els @cmla . ens-cachan . fr 
Centre de mathematiques et de leurs applications 
Ecole normale superieure de Cachan & CNRS, UMR 8536 
61, avenue du President Wilson 
94235 Cachan cedex, France 



In this paper, we study a simple one-dimensional model of reaction- 
diffusion with bistable non-linearity and in heterogeneous media. The 
bistable term accounts for the so-called Allee effect, and the hetero- 
geneity in the media is localized. We recall the definition of a transi- 
tion wave used in similar situations by well-known authors, and pro- 
pose an alternative definition for discussion. We call it generalized 
traveling wave. As a consequence, we give new results of pinning in 
such media, due to Allee effect. 
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1 Introduction 



1 . 1 Framework 

We will focus in this paper on the following reaction-diffusion equation: 

d t u = Au + r(x) f(u) 
u(-,0) =u 

It models a wide range of phenomena in biology, which will be our main 
concern in applications. While the homogeneous version of this equation 
(r(x) = r ) has been well studied since the pioneering works of Fisher [ 
and Kolmogorov and al. [3], the heterogeneous ones are an actual field of 
research. We assume here that heterogeneity occurs only in the reactivity 
of the equation, through the r(x) function, while the diffusivity is constant 
(d = 1 can be achieved by a convenient rescaling). 

The most part of the paper is devoted to the so-called bistable reaction 
term: 

f(u) = u(l-u)(u-a) (2) 

with < a < 1. It is even convenient to suppose < a < |, as we will see 
below. Whenever possible, we will discuss the extension of our results to the 
monostable case (Fisher-KPP) 

f(u) = u{l - u) (3) 

Our interest here is in equations that are local perturbations of the ho- 
mogeneous one: we suppose then that r(x) = r + r 1 (x), with r > 0, r 1 
compactly supported, continuous, and so that r(x) remains non- negative: 

I In I |oo < r . 

While it is very common to see applications to biology make use of a 
monostable reaction term, let us recall that in a number of cases it is not a 
reasonable assumption near 0. In our mind, the reaction term describes the 
growth of a biological quantity (be it a population density, a concentration, a 
proportion of a population sharing a given attribute, ■ ■ ■ )• A logistic growth 
is the first improvement over an exponential one (f(u) = ru), taking account 
of the existence of the habitat or environment (with capacities, resources, 
limits). Under such an hypothesis, the growth is no more proportional to u, 
and saturates when the population reaches the so-called carrying capacity of 
the habitat K > (f(u) = ru(l - f )). 
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The next improvement concerning the reaction term is related to the 
Allee effect. The Allee effect is the correlation between individual fitness and 
global density. At low densities, the survival of an individual can be made 
more complicated by increased predation or difficulties in mating. In this 
kind of situation, it it not reasonable to have / logistic at low densities. We 
can differentiate between weak Allee effect (for which /'(0) = is a good 
approximation), and Strong Allee effect (/'(0) < and thus / < near 
zero). The bistable term models then a population growth with a strong 
Allee effect. Such a case is also called the "heterozygote inferior" in a genetic 
context [ ]. 

The reactivity r takes account of the intensity of the reaction, here inten- 
sity of the growth (population case), of the selection (genetic case). 



1.2 Classical theory 

We recall here quickly the known results for equation (I) when r is constant 
equal to 1 and / bistable: 

d t u = Au + u(l — u)(u — a) 
u(-,0) = u 

We are interested here in classical traveling waves of this equation: u(x, t) = 
U(x — ct). For clarity, we choose to investigate waves going to the right, with 
invasion of the state by the state 1. Such a U is solution of 

U" + cU' + U(1-U)(U -a) U(—oo) = 1 U(+oo) = (5) 

It is known [4, 5] that U is unique up to translation, and c is related to a: 

U(z) = — c=\-(l-2a) (6) 

v ; 1 + exp (z) V 2 v ; v ; 

We see that the speed c is positive as long as a < |. This is related to the 
sign of the integral of f(u) on [0, 1]. For a constant positive r we have 

U{z) = — \ 7^ c=J\(\-2a) (7) 
I + exp [Zy/l) V 2 

A comment that can be done here is the r-dependence on the speed and 
profiles of the wave. Formulas in (7) show that bigger r leads to quicker and 
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steeper wave. Thus we can say that when r is big, the media is favorable for 
the invasion by the state 1 and the invasion is sudden. When r is low, the 
media is less favorable, and the invasion is progressive. 

2 Heterogeneous media: transition waves vs. 
generalized traveling waves 

2.1 Transition waves 

More general media can exhibit range more or less favorable for propagation. 
As explained before, we study here media that are local perturbation of the 
homogeneous case. This allows the relaxation of the initial condition to the 
traveling wave, due to convergence results (see for example [5]). Then, from 
an initial condition far enough from zero, the solution is near the classical 
traveling wave before it reaches the zone where r varies. 

We expect that in the bistable case, stability of the leading edge of the 
wave will lead to similar behaviors when r varies. This has been reinforced 
by numerical simulations made in [6]. 

To investigate such situations, one can use the notion of transition fronts 
defined for example in [7] 1 . A transition front for equation (1) is a global in 
time solution that connects the two steady states: 

u(—oo,t) = 1 u(+oo,i) = 

for each t, and whose e-width is uniformly bounded in time for all e. The 
£-width is defined as follow: diam{x G M.\u(x,t) > e and 1 — u(x,t) > e} 
This has been done in [8, 9]. 

2.2 Generalized traveling waves 

We propose here another approach, already developed in [6], though we did 
not state there an appropriate definition. The main idea is to follow the front 
at each time (each time step for the numerics). We then get another type of 
result: an approximation result instead of a convergence one. 

1 It seemed to us during our bibliographic research that these authors previously used 
the term of generalized traveling wave, before they coined the term transition front (or 
transition wave). They seem now to stick to the latter. 
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Our work is backed upon an extensive numerical study, which gave birth 
to a conjecture. This conjecture reads: every reasonable initial datum will 
evolve following the scheme 

where U is the solution of the homogeneous equation (4), x the center (or 
position) of the wave, w its width. This width is defined in the same way 
as above: diam{x G M.\u(x,t) > 5 and 1 — u(x,t) > 5} for an appropriate 
S chosen once for all. The profile U makes it unnecessary to control the 
width for any 5, for widths for different thresholds are related with a fixed 
coefficient. 

Conjecture 1: 

Let u be the solution of the Cauchy problem (1) with reaction term (2). 
Let U be the wave profile in homogeneous medium defined in (6). 

Then there exist Xq and w such as u(x,t) — U remains small in 

L°°-norm provided as r varies few enough: 

We > 3r] > such as Vt > 3(x (t),w(t)) 
INU < V => IK,*) - u C~J^ ) Hoc < £ 

Definition 1: 

Such a solution of (l)-(2) is called a generalized traveling wave 2 . 

Full proof of this conjecture will be left for a future paper. We want to 
emphasize here that this definition covers cases of defects narrow or wide, 
symmetric or antisymmetric as defined in [6] . The big defects though are not 
covered, but they behave in a similar manner. 

If this conjecture holds, and we strongly believe it to do so, it is an 
interesting result because we have ways to compute x and w. This is dealt 
with in the next section. 

Note also that this conjecture does not hold for the Fisher-KPP case 
due to the instability of the state 0. Instead, secondary outbreaks occur in 
advance of the wave [6]. 

2 We use this term because we believe it is left free by the authors previously mentioned. 
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3 Velocity and width 
3.1 Moments 

Given that U(z) = tt^ — , we have U'(z) = - t— 

v / 1+expz' V / (l+cxpjz)- 2 

decreasing at infinity and its moments are well defined. The three first ones 



cxp z 



r. Thus U' is fastly 



are: 



' r+oo 

oiq = U'(z) dz = — 1 

J — oo 



Oil 



>3C 

+oo 



zU'(z) dz = 



(8) 



■oo 
+oo 



If we assume very roughly that we have 

x — Xo(t) 



y 



u(x, t) = U 



w(t) 



then we can differentiate u: 

d x u(x,t) 



1 , ( x- x (t) 



w(t) V w (t) 



Making, for a given t the change of variable z = x in the three integrals 



of 



yields 



+oo 



Xq 



W = — - 



■oo 
+oo 



71- 



d x u(x, t) dx 
xd x u(x, t) dx 

x 2 d x u(x, t) dx + x, 



(9) 



So that we have formulas for x and w. This first way to obtain them is the 
least interesting, because it requires to compute (numerically) the solution 
u, then its partial derivative, and some additional integrals. 
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3.2 Least square fit 



The second way to obtain the position and width of the generalized traveling 
wave has been developed in [6]. It consists in fitting (at each time) the 
solution u with an appropriate profile. We have consistently checked that 
this fitting procedure can be done with a numerical routine providing a very 
small error (in L°°-norm). This evaluation of xq and w is excellent, and takes 
less time than the one above. 

3.3 Balance laws 

A much better solution is given by the following system (10). This system 
comes from two simple balance laws from the original partial differential 
equation (1), and has been derived in [6]. It is a basic system of two ordinary 
differential equations that yield the evolution of the position and width: 



By this way, evaluation of the position and width is simpler and quicker. 

4 Applications 

We derive in this section some easy consequences of the equations (10). These 
equations have the great benefit to link the position and width of the gener- 
alized traveling wave to the heterogeneity r. They show that the qualitative 
behavior of the front around an obstacle is not trivial (see [6] for details). 
We have notably a slowing down of the front before it reaches any positive 
defect. We add here two more pinning criteria to the one previously given in 
this article. 

4.1 Pinning of the front by an antisymmetric defect 

Let us consider the case of a narrow antisymmetric defect, modeled by an 
Heaviside step function: 




(10) 



r(x) = r + ri H(x) 
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where H(x) = if x is negative or superior to a given constant R big enough, 
and H(x) = 1 between and R. We deal with the case T\ > 0. 

1. Equations for the critical width 

The system of equations (10) can be rewritten for such a defect r(x), 
and when the front stops we have: 



we get: 



r Q w 



1 - 2a 



1 

3w 



w 



or: 



-r - 



r Q w 



w 



+riw 



+riw 



0, w 





f{U(z)) dz 

(l-2U(z))f(U(z)) dz 



1 - 2a 



To — 

3w 6 



+riwSi 
+r 1 wS 2 



w 

-Xo 

w 











(11a) 
(lib) 



We can first simplify the first equation by w, then we note that the left 
part of the left hand side of (11a) is always positive. In the right hand 
side, we supposed r% > 0. To balance the equation, it is necessary that 
Si (^f-) contributes negatively. This is indeed possible given the plot 
of Si(y) = ^i+Jj^P '■ Si has a negative minimum (see figure 1). 




Figure 1: Plot of Si. We chose a = 0.25. 
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The front stops if and only if: 



2 

For the critical case, we have equality: 

r ° ("~ 2~~) = -ricmin (' S ' 1 ) ( 12 ) 

This shows the existence of a threshold over which pinning occurs, and 
under which the wave goes on propagating. Note that this pinning is 
counter-intuitive for r\ > 0, given the considerations before given on 
the r-dependence of the speed. 

2. Resolution of the equations 



The minimum of S\ and its argument z m i n can be obtained because 

(1— a— a< 
(l+e»)3 



S[(y) = -R x {y) where R x {y) = f(U(y)) = "^£f\ We have 



Ri(y) = o 1 - a - ae y = 0, so 

Zmin 111 I 

\ a 
then 

. ,<n or ^ l-2a(l + ^) a 2 
min(Si) = S'll^mm) = + = -y 

This gives the critical height thanks to (12) : 

r (l - 2a) 



ric 



a 2 



We find next the width w at pinning thanks to (13b), with the substi- 
tution of = ^ a by In (^^) and of r\ by the obtained quantity. After a 
little algebra, we get : 



1 

w 



r v /l + (l-2a)(3-2a) 
and finally : 



r v /l + (l-2a)(3-2a) 
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We have evaluated this result using the following numerical application 
with a = 0, 3 and r = 1: 



' 40 

r lc = — « 4.44 
9 

'50 
49 




,x = -W — In - « -0.86 



which compares very well to the values coming from the simulations (error 
of 3%). 

4.2 Pinning of the front by a no-reaction zone 

Let us consider the case of a no-reaction zone, modelized by the following 
function: 

{r if x < 
if < x < X\ 
r if x > X\. 

We ask the question of the minimal width required to stop the front. We call 
this minimal width x± c , so that we have the following dichotomy: 

• if X\ < Xic, then the front is not stopped by the obstacle. We will say 
that the front goes through the obstacle, it goes on propagating after 
the zone. 

• if x\ > Xic, then the front is stopped, pinned in the no-reaction zone. 
The wave is not moving any longer, neither in one direction nor in the 
other: its speed of spread decreases to zero. When it reaches the no- 
reaction zone, the front slows down and changes shape. At the limit, 
the front tends to an aifine profile (solution of the heat equation with 
Dirichlet boundary conditions). 

We solve this problem in a straight-forward manner: 

1. Equations for the critical width 

The system of equations (10) can be rewritten for such a defect r(x), 
and when the front stops we have: 

xq = 0, w = 
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Then we get: 



w r 
1 

it- + w r 

or stated otherwise: 

w r (ti 
1 



f{U(z)) dz + 



+00 



x-i^ — xrt 



f(U{z)) dz 



[l-2U(z))f(U(z)) dz 



X-tr—Xn 



-Xq 



3w 



+ w r [T 2 



-Xq 
W 



X\ c — Xq 



W 



X\ c — Xq 



■w 



[l-2U(z))f(U(z)) dz 
(13a) 

= 0. (13b) 



The first equation can be divided by w r (positive quantity), and the 
wave stops when both terms in the left hand side of (13a) balance each 
other. The integrals can be calculated, we have: 



Tx(y) 



e y(2 + e y - 2a(l + e»)) 

2(1 + e y) 2 



It is obvious that Ti(y) = ^^2n^^^^ ls a l wa Y s positive for all 
a G [0, |]. We show in figure 2 an example of such a plot. 




Figure 2: Plot of T x . We chose a = 0.25. 

So that the equation can be verified, it is necessary that Si ( Xlc ~ x ° ) be 
negative. This is indeed possible given the profile of Si (see figure 1). 
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S% has a strictly negative minimum, but for widths below the critical 
width, this negative part is not enough to stop the front. For the critical 
width x\c, we have exactly 



w 



T 1 [^)=-mmS 1 (14) 



2. Values of and Xlc x ° when the wave stops 

WW 1 

We have seen previously that the minimum of S% and its argument z ri 
can be easily computed: 



In 



1 



a 



min(S'i) = S x (z m i n ) = 



2(1 + 



We have then 



xic - %o _ , (I - a 

Zmin m 

W 



and we get the expression, as a function of a and b = \fl — 2a, of the 
ratio of both collective coordinates when the front stops for the critical 
width: 

Xq ( a + a 2 — 1 + b s 

— m 



w \ 2a + a 2 — 1 

thanks to equation (14), because — min(S'i) has a unique antecedent 
byTx. 

3. Resolution of the equations 

We replace the quantities = ^ Q - and Xlc ~ x ° by their values in equation 
(13b), we can then solve and we obtain the width w c of the front when 
it stops: 

(b - aY 

w c - 



JWq V -1 + 3a - 8a 3 + 9a 4 + 7a 4 6 -2ab + b- 3a 5 - 2a 3 6' 
and the critical width well. The pinning criteria writes: 

x\ ^ x lc = w c ( In ( - — -j-lnf-^- 1 ^ — 1 - b ) ) (15) 



2a + a 2 - 1 
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It is interesting to note that the critical width depends on r via the inverse 
of its square root. We can check on figure 3 that the right hand side of (15) 
is a decreasing function of a on [0, |], with infinite limit at and null at |. 

30 
25 
20 
15 
10 

5 



0.1 0.2 0.3 0.4 0.5 

Figure 3: Plot of ^ lM ^^ 2ab+b .^_ 2cfib (in (^) - In 

This is a proof of consistency of the standardized model, as we expect the 
front to be more easily pinned as a approach the value \. 

The values obtained via this criteria are in very good agreement with 
numerical simulations of the initial partial differential equation. For r = 1, 
our criteria (15) gives X\ c ~ ~ 7.18 (we recall that we fixed a = 0.3). 
Numerically, we find that the front stops from x\ c ~ 6.90, which yields an 
error inferior to 4%. 

5 Conclusion 

In this paper, we proposed a notion that generalizes the classical traveling 
wave in homogeneous media to heterogeneous ones. It uses two moments of 
the derivative of the profile U, which has the appropriate behavior at infinity 
for doing so. 

We have quickly shown in this paper the power of this new tool. Once 
written the evolution equations of the position and width, it becomes possible 
to tackle down difficult problems while answering questions that are of crucial 
interest to biologists. The idea of following the front at each time is fruitful for 
people interested in applications, and complements the convergence results 
obtained in a traditional manner. 
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